function [ yi ] = Interp( x, Y, xi, method, varargin)
%Summary:  用method指定的插值方法做一维插值，输出点列xi的数值
%变量说明：
%   x 横坐标数组,要求互不相同。
%   Y 纵坐标数组
%   xi 求此点列的数值
%   method，插值方法
%       method 包括：
%           1，'linear'      线性插值
%           2，'newton'      牛顿插值
%           3，'lagrange'    Lagrange插值
%           4，'hermite'     Hermite插值
%           5，'spline'      三次样条插值
%    输出变量：yi表示xi点列对应的数值。
    if length(x)~=length(Y)
        error('输入数据错误！')
    end
    n = length(x);
    ln = length(xi);
    
    %下面给x排序。
    [x,IX] = sort(x);
    for i = 1:n-1
        if x(i)==x(i+1)
            error('数据错误，x中有相同的值');
        end
    end
    Y = Y(IX);
    
    yi = zeros( 1, ln);
    
    %线性插值：
    if strcmp('linear',method) == 1
        r = 1;
        
        %先计算各线段斜率。
        k = zeros(n,1);
        for i = 1 : length(x)-1
            k(i) = (Y(i+1)-Y(i))/(x(i+1)-x(i));
        end
        
        for i =1:ln
            while x(r)<xi(i)
                 r = r + 1;
            end
            if r == 1
                yi(i) = Y(1);
            else
                yi(i) = Y(r-1)+k(r-1)*(xi(i)-x(r-1));
            end
        end
    end
    
    %Newton插值 
    
    %递归求差商。     
    function [y] = Dif(i,j)
        if isnan(dif(i,j)) && j > i
            dif( i, j)=( Dif( i, j-1) - Dif( i+1, j))/( x( i) - x( j));
        end
            y = dif( i, j);
    end

    if strcmp(method, 'newton')==1
        %先求出差商dif(i,j)=f[x(i),x(i+1),...,x(j)]
        dif = nan(n,n);
        for i = 1:n
            dif(i,i) = Y(i);
        end
        
        for i = 1 : ln
            yi(i) = Y(1);
            p = 1;
            for j = 2 : n
                p = p*(xi(i)-x(j-1));
                yi(i) = yi(i) + Dif(1,j)*p;
            end
        end
    end
    
    if strcmp(method,'lagrange') == 1
        %先求各基函数分母
        fenmu = ones(1,n);
        for i =1:n
            for j = 1:n
                if j ~= i
                fenmu(i) = fenmu(i)*(x(i)-x(j));
                end
            end
        end
        
        for i = 1:ln
            T = 1;
            for j = 1: n
                T = T * (xi(i)-x(j));
                if T == 0
                    yi(i) = Y(j);
                    break
                end
            end
            
            if T == 0
                continue
            end
            for j = 1: n
                p = T/(xi(i)-x(j));
                p = Y(j)*p/fenmu(j);
                yi(i) = yi(i) + p;
            end
            
        end            
    end
    %hermite插值
    
    if strcmp(method, 'hermite')==1
        h = x(2:n)-x(1:n-1);
        if nargin < 5
            error('请指定一阶导数')
        end
        Y1=varargin{1};
        r=1;
        for i = 1:ln
            while x(r)<xi(i)
                r = r+1;
                if r>n
                    yi(i)=Y(n)
                    continue
                end
            end
            if r == 1
                yi(i) =Y(1);
                continue
            end
            hr = xi(i) - x(r);
            hl = xi(i) - x(r-1);
            pr = ( hr/h(r-1))^2;
            pl = ( hl/h(r-1))^2;
            yi(i) = Y(r-1) * pr*(1+2*hl/h(r-1));
            yi(i) = yi(i) + Y(r)*pl*(1-2*hr/h(r-1));
            yi(i) = yi(i) + Y1(r-1)*hl*pr;
            yi(i) = yi(i) + Y1(r)*hr*pl;
        end        
    end
    
    %自然三次样条
    if strcmp(method, 'spline')==1
        
        h = x(2:n)-x(1:n-1);
        b = 6*(Y(2:n)-Y(1:n-1))./h;
        u = 2*(h(1:n-2)+h(2:n-1));
        v = b(2:n-1)-b(1:n-2);v=v';
        A = zeros(n-2);
        for i = 1:n-2
            A(i,i) = u(i);
        end
        for i = 1:n-3
            A(i,i+1) = h(i+1);
            A(i+1,i)=A(i,i+1);
        end
        z = zeros(1,n+1);
        z(2:n-1) = (A\v)';
        
        h2 = h.^2;
        c = (Y(2:n) - z(2:n).*h2/6)./h;
        d = (Y(1:n-1) - z(1:n-1).*h2/6)./h;
        
        for i = 1:ln
            r = 1;
            while x(r)<xi(i)
                r = r+1;
                if r>n
                    yi(i)=Y(n)
                    continue
                end
            end
            if r == 1
                yi(i) = Y(1);
                continue
            end
            l = r-1;
            hr = xi(i)-x(r);
            hl = xi(i)-x(l);
            yi(i) = c(l)*hl-d(l)*hr;
            yi(i) = yi(i) - z(l)*hr^3/(6*h(l));
            yi(i) = yi(i) + z(r)*hl^3/(6*h(l));
        end
        
    end
end

